Molecular line radiative transfer in protoplanetary disks: 
Monte Carlo simulations versus approximate methods 
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O ■ ABSTRACT 

^ i We analyze the line radiative transfer in protoplanetary disks using several 

I approximate methods and a well-tested Accelerated Monte Carlo code. A low- 

mass flaring disk model with uniform as well as stratifled molecular abundances 
is adopted. Radiative transfer in low and high rotational lines of CO, C^^O, 
HCO+, DCO'*', HCN, CS, and H2CO is simulated. The corresponding excitation 
temperatures, synthetic spectra, and channel maps are derived and compared 
to the results of the Monte Carlo calculations. A simple scheme that describes 
the conditions of the line excitation for a chosen molecular transition is elabo- 
rated. We flnd that the simple LTE approach can safely be applied for the low 
molecular transitions only, while it significantly overestimates the intensities of 
the upper lines. In contrast, the Full Escape Probability (FEP) approximation 
can safely be used for the upper transitions (Jup > 3) but it is not appropriate 
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for the lowest transitions because of the maser effect. In general, the molecular 
lines in protoplanetary disks are partly subthermally excited and require more 
sophisticated approximate line radiative transfer methods. We analyze a number 
of approximate methods, namely, LVG, VEP (Vertical Escape Probability) and 
VOR (Vertical One Ray) and discuss their algorithms in detail. In addition, two 
modifications to the canonical Monte Carlo algorithm that allow a significant 
speed up of the line radiative transfer modeling in rotating configurations by a 
factor of 10-50 are described. 

Subject headings: astrochemistry - line: formation; profiles - radiative transfer 
- planetary systems: protoplanetary disks - stars: formation - radio-lines: stars 



1. Introduction 



The study of the formation and evolution of protoplanetary disks is an important topic 
for our understanding of the formation of extrasolar planets and the planets in our Solar sys- 
tem. The gas content of the disk that represents most of its mass and has an imprint of earlier 
evolutionary history is best studied in molecular lines. Since 1997, a number of rotational 
transitions from simple molecules like CO, CS, CN, HCN, HCO+, N2H+, H2CO, C2H, and 
some of their isotopes have been firmly detected in several nearby protoplane t ary disks with 



single-dish telescopes and millimeter interferometers fe.g.. iDutrev et al 



Guilloteau et al. 



19971: Ivan Zadelhoff et all I2OOII : lAikawa et all l2003l : lOi et al 



200 



i 



199^: iKastner et al.l 



2004 



Thi et al. 



2004 



20061 ). The wealth of information about disk physical structure and chem- 



ical composition can only be acquired with multi-molecule, multi-line interferom etric ob 



servations at (sub-) mi llimeter wavelengths and at sub-arcsecond resolution (e.g., lOi et al. 



20061 : iPietu et al.l 120071 ). In 2005, we started a comprehensive observational program at the 
IRAM 30-m antenna and the Plateau de Bure interferometer, aimed at studying protoplan- 
etary disks with several key molecular tracers of disk temperature, density, kir iematics, and 
ioniz ation degree ("Chemistry in Disks" (CID) project; for the first results see iDutrey et al. 
2007h . 



However, in order to extract all the relevant information from molecular lines and 
thermal dust continuum one has to cope with the inverse (line) radiative transfer (LRT) 
problem, which can hardly be solved for complex geometrical configurations. Instead, one 
solves the direct LRT problem: models with varied parameters are compared to the ob- 
serva tional data until the best-f i t is found by using some iterative m inimization technique 



[e.g. iGuilloteau fc Dutreyl Il998l : iDartois. Dutrey. fc Guilloteaul l2003l ). Naturally, this task 
requires enormous amounts of computing power that severely restricts the choice of applica- 
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ble computational methods. This imphes that often fast but only approximate algorithms 
can be used. Thus it is crucial to investigate the applicability of such approximate methods 
for the LRT modeling of protoplanetary disks. 

Various one- and multi-dimensional exact and approximate LRT methods have been 
developed for modeling of stellar atmospheres, stellar winds, interstellar clouds, and star- 
forming regions (for a good overview see, e.g., iPeraiah 1 12004| ). Among the simplest LRT 
concepts are the LTE (Local Thermodynamical Equilibrium) and LV G (Large Velocity Gra- 
dient s) approaches in which the LRT problem is treated locally (see ISobolevlll960l : iMihalas 
19781 ). The LVG method is just one among the sub-class of "Escape Probabihty" approaches 
that allow the straightforward calculation of the escape probability. These simplified LRT 
methods have been applied to constrain the basic ph ysical parameters o f dense molecular 
cores using single- dish mole cular line data (see, e.g. 
198ll : iGoldsmith et al.lEgSsl ). 



de Jong et al.lll975l : iGuilloteau et al. 



Later more sophi sticated, non-local LRT approaches have been developed, starting from 
the ID methods (e.g.. lLucaslll974l : lElitzur fc Asensio Ramosll2006l ) to the multi-dimensional 
(Acceler ated) Monte Carlo, Local Linearization, and Multi-Zone Escape Probability meth- 



ods, e . g. iBerned (119791) 



f 2000 ): bssenkoDf et al 



2001) 



Juvelal ( 19971) :lDullemond &: TuroUal (120001) :lHogerheiide &: van der Tak 



Rawlings fc Yated (1200111 : ISch5ier fc OlofssonI fl200ll\: 



Keto et al. 



hooi ): IPoelman fc SpaanJ ( I2OO5 I) (for a recent review see Ivan Zadelhoff et al.ll2002l ). Being 



more accurate, modern LRT methods are routine ly applied to study the kinematics and 



Pavlyuchenkov et al.l 20061: Tafalla et al. 



chern ical structure of prestellar cores in detail (e.g. 
20061). but their us e for protoplanetary disks is still the exception (e.g. ISemenov et al.ll2005l : 
Brinch et al1l2007h . 



The global appearance of the physical structure of protoplanetary disks is commonly 
treated in the framework of steady-state passive or active flaring disk models with or without 
an inner rim, which seem to be successful in ex plaining many observational facts (e.g.. 



D'Alessio et al.lll998l . 1200 ll : iDuUemond et al.ll200ll ). The modeled structure can further be 



used to unravel the chemical evolution of the disk with a sophisticated gas-grain chemical 

e^ 



Willacv et al. 


19981 


^^^^ f^^^^-^^^^^ 
lAikawa et al. 


2OO2I: 


Gail 


2OO2I; 


2006|; 


Willacv et al. 


2006|; 


Tscharnuter & Gail 


2OO7I) 



3; lllgner et all booi ISemenov et"all 12004 



2007h . Finally, the disk physical structure and 



chemical composition together with appropriate collision rates provide the necessary inputs 
for the LRT simulations. 

Since the iterative procedure to search for the best fit to the observed spectra can easily 
become computationally prohibitive when the full LRT has to be considered, it is of great 
importance to identify those fast approximate LRT methods that give sufficiently accurate 
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results (with respect to the observational hmitations) and to investigate the hmits of their 
apphcabihty. The analysis is further complicated when a realistic disk chemical structure 
is taken into account. To the best of our knowledge, a detailed benchmarking of various 
approximate LRT approaches for p rotoplanetary disks has no t yet been attempted in the 



literature apart from the work by Ivan Zadelhoff et al.l (120011 ) and therefore constitutes a 
major goal of the present paper. 

In this paper, we compare a full 2D non-LTE Monte Carlo code with several approximate 
line radiative transfer methods by considering a few representative radiative transfer cases. 
These problems are thought to cover all possible situations encountered in disks, namely, the 
excitation of optically thin/thick emission lines in low/high density regions. Consequently, 
we consider low and high transitions of key molecules in disks tracing temperature (CO and 
isotopologues), column density (^'^CO or C^^O), density (CS and H2CO), deuterium fraction- 
ation (DCO^), ionization degree (HCO"*"), and photochemistry (HCN). For this exploratory 
study we adopt an irradiated flaring model of a young, low-mass T Tauri protoplanetary 
disk, and use both uniform molecular abundances and those from an advanced chemical 
model (Section ED. Our 2D Monte Carlo algorithm is based on the ID code "RATRAN" 



( iHogerheijde fc van der TakI l2000l ) with two acceleration modifications for disks (see Sec- 
tion [3]). The adopted approximate methods include the classical LTE approach, the Full 
Escape Probability (FEP) and Vertical Escape Probability (VEP) methods, LVG, and a 
non-local ID-method VOR (Vertical One Ray method). The approximate methods are all 
described in Section m Using the input disk model and the above LRT methods, in Section [S] 
we calculate the corresponding excitation temperatures, synthetic spectra, and channel maps 
and compare these data to the results obtained with the Monte Carlo code. The discussion 
and final conclusions follow. 



2. Disk physical and chemical structure 



We consider a typical protoplanetary disk surrounding a T Tauri-like star in the nearby 
Taurus-Auriga molecular complex at 140 pc. A flaring stead y-state model represen tative of a 



Class II disk with vertical temperature gradient is adopted (ID'Alessio et al. 



19991 ). The disk 



has inner and outer radii of 0.037 and 800 AU, respectively, a mass of 0.07 M0, and a mass ac 
cretion rate of 10~^ Mq yr~^. The central star has a radius of 2.64 Rq, a mass of 0.7 Mf^, an ef- 



fective temperature of 4 000 K, and an age of about 5 Myr (see Table 2 in lDutrey et al.ll2007l ). 
We assume that the disk has a constant m icroturbulent velocity of 0.1 kms~^ a nd a purely 
Kepl erian rotation profile, V{r) oc r~^/^ (IDartois. Dutrey. k. Guilloteaul l2003l : IPietu et al. 



20071 ). The disk structure is shown in Fig. [T]and the parameters are summarized in Table [TJ 
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In the present study we focus on the molecules that have been detected in protoplanetary 
disks and are readily used as probes of the disk physical and chemical structure, high-energy 
radiation, and deuterium fractionation: CO, i^CO/C^^O, HCO+, DCO+, HCN, CS, and 
H2CO. Corresponding to the different excitation conditions, the various rotational lines of 
these molecules are thought to be generated in distinct disk regions ranging from the cold 
midplane and warm molecular layer to the hot disk surface layer. 

Normally one has to use a sophisticated gas-grain chemical network with surface re- 



the disk (see 


van Dishoeck & Bla 






1998 


Semenov et al. 


2005: 


Willacv et al. 


2006). 



1998; Aikawa fc Herbst 1999; Markwick et al. 2002 



are reached in the disk intermediate layer (at ~ 0.5 — 1 pressure scale heights) that is suffi- 



Aikawa et al. 


2002; 


Dutrev et al. 


2007) 



models, from a qualitative point of view the abundance distributions of many molecules can 
be represented by a single layer with particular concen tration, thickness, and location in the 



disk (see, however. ISemenov et al.ll2006l : lAikawal 120071 ). 



However, for the first simple analysis of the observational data it is often sufficient to 
assume that molecular abundances are constant (with respect to hydrogen). Moreover, the 
use of uniform molecular abundances can help to disentangle most important line radiative 
transfer effects from strong chemical gradients in the disks. 

We adopt both approaches. The values of the uniform molecular concentrations and 
parameters of the molec ular layers are o btaine d by using the gas-grain chemical model 
with surface reactions of ISemenov et al.l (120051). In essence, this model is based on the 



UMI ST95 database of gas-phase reacti ons (IMillar et al.lll997h. a set o f surface reactions 
from Hasegawa. Herbst. fc Leunej Jl992h and iHasegawa fc HerbstI Jl993h on 0.1/ x m sili cate 
grains, and includes a small deuterium network from lBergin. Neufeld. &: Melnickl ( 119991 ). 



The location of each molecular layer in the disk is constrained by the upper and lower 
boundaries that can be approximated by a power law: z(r) / HJr) axr^ (where r is the ra- 



dius in AU and Hp{r) is the pressure scale height as defined in lDartois. Dutrey. fc Guilloteau 



2OO3I ). We assume that the molecular abundances are constant across a molecular layer and 
utilize the maximum values in the uniform disk model (Fig. [1]). Note that the CS abun- 
dances peak toward the inner disk region and cannot be accurately represented by such a 
simple one-zone model. The parameters of the uniform and layered chemical structures, i.e. 
parameters flmm, &mm and Umax) bmax describing the upper and lower boundaries as defined 
above, are presented in Table [2l 
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Finally, the collision rates for the line radiative transfer mo deling; are take r i from the 



pubhcly available "Leiden Atomic and Molecular Database" 0, see Schoier et al. ( 2005 ). 



3. Accelerated Monte Carlo method 

In this section we describe our implementation of the Accelerated Monte Carlo technique 
in the 2D non-LTE line radiative transfer code ART, which is part of the software package 
"URAN(IA)" aimed at modeling of various objects with different exact/approximate LRT 
methods. 



3.1. Basic algorithm 

The general idea of all Monte Carlo radiative transfer methods is to follow the propaga- 
tion of randomly ejected photons or photon packages through the medium. These codes are 
based on a straightforward formalism that is easy to implement numerically and are capable 



complex multidimensional ob 


ects with arbitrar- 


Juvela 


1997 




Wolf et al. 


1999 




Pinte et al. 


2006; 



NicoUini fc Alcolea II2006I ). However, in order to lower random noise and thus retain reliable 
accuracy of the results, Monte Carlo simulations typically require enormous amounts of com- 
puting time, especially for highly optically thick objects. In certain cases the numerically 
most demanding part of Monte Carlo codes (ray tracing) can be significantly optimized by 
special iterative schemes for the acceleration of convergence. 

The algorithm of our 2 non -LTE line radiative transfer code ART is fully described in 
Pavlyuchenkov fc Shustovl (120041 ) and only briefly summarized here. In essence, this Monte 
Carlo code is based on the accelerated scheme originall y prop osed and implemented in the 
LRT code "RA TRAN" by iHogerheiide fc van der Ta^ J2000h and in the modified Monte 
Carlo method of Voronkov ( 1999h. The rnain equations of the line rad iative transfer can be 
found elsewhere (e.g.. Ivan Zadelhoff et al.ll2002l : iRohlfs fc Wilsonll2004l ) and are not repeated 
here. 

First, the computational domain is split into cells. All physical quantities are assumed 
to be constant within each disk cell i, except for the velocity that can vary. In the first 
iteration, initial level populations and a set of photon random directions n{i) through the 
grid are defined (this is the Monte Carlo part of the calculations). The same set of photon 



^http : //www ■ strw . leidenuniv . nl/$\siin$mold.ata 
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random directions is used in all subsequent iterations for accurate ray tracing, producing a 
solution free of random noise but with possible undersampling of directions and velocities. 

Using these values, the line intensities luii) are computed for each disk cell i by explicit 
direct integration of the radiative transfer equation with the long characteristic method along 
the predefined set of photon rays n{i), starting with the CMB field at the edge of the disk (or 
interstellar radiation field, or anisotropic stellar radiation). The integration path is split into 
sub-intervals with a step that is smaller than the cell size. Ideally, one would need to define 
this integration step based on local gradients of the physical conditions in the disk and the 
optical depth of the line. However, this is not easily achievable in practice, so an empirical 
value has to be used. We adopt a step size of 1/10 of the cell size in our computations. 

Next, the resulting mean line intensity J is calculated in each disk cell by spatial av- 
eraging of the derived intensities over directions n{i) and by averaging over the line 
profile. The computed mean intensities are utilized in the next iteration step to refine the 
level populations by solving the balance equations in all disk cells. 

The convergence of the integration procedure for optically thick lines can be substan- 
tially accelerated by ad ditional internal sub-it erations similar to the Accelerated Lambda 



Iterations (ALI) scheme (lAuer &: Mihalaslll969l ). The ALI method employs the fact that the 



calculated mean line intensity at each particular grid cell can be separated into an internal 
component generated in the cell itself and an external contribution coming from other grid 
cells. As a result, additional sub- iterations are applied to bring into an agreement the in- 
ternal mean line intensity and the corresponding level populations in every grid cell. When 
these local sub-iterations converge, the global iterations over the entire grid are continued 
until the level populations in each cell do not differ by more than 0.1%. Note that this 
scheme may not be fully appropriate in cases of particularly thick lines with r > 100 when 



such a simple convergence criterion fails, e.g. for water lines (Ivan der Tak et al.ll2005l ). 



After the final molecular level populations are obtained, the maximum random error 
associated with the Monte Carlo computations is estimated. For that all iterations are 
repeated again with another set of preselected random photon rays. We choose the number 
of random photon paths such that the relative error in the computed level populations 
is smaller than 1%. Although the code is currently restricted to axially symmetric (2D) 
radiative transfer problems in polar coordinates, it can be extended to problems with any 
number of dimensions and for any coordinate system. 

Finally, these resulting level populations are used in the 3D ray tracing part of the 
"URAN(IA)" package to produce synthetic molecular line spectra, channel maps, integrated 
line intensity maps, velocity maps etc., with an arbitrary number of velocity channels, and 
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for any disk viewing geometry. Beam convolution can be applied if needed. 

The "URAN(IA)" code was carefully benchmarked against other LRT codes by solving 
several one- and two-dimensional problems. In the one- dimensional regime, the computed 
molecular level populations and excitation temperatures are fully consistent with all tests 
proposed during the confe rence on molecular radiative transfeiB (Leiden University, 1999, see 



van Zadelhoff et al.ll2002l ). Several additional tests of our code were performed in 2D mode. 



using a few simple problems with analytically predictable solutions. In addition, together 
with M. Hogerheijde - the developer of the two-dimensional LRT code "RATRAN" - we 
successfully compared the results for both codes for Keplerian disks. 



3.2. Further acceleration of the Monte Carlo algorithm 

3.2.1. Concept of interaction areas 

The method of long characteristics is one of the most reliable algorithms to calculate the 
mean intensity and thus to o btain an accurate solution to the radiative transfer problem (e.g.. 



DuUemond fc Turollall200d ). However, a direct realization of this method is computationally 
expensive for multi-dimensional radiative transfer simulations. For example, to compute 
molecular level populations for a 2D disk model having 150 vertical and 50 radial grid points 
and 1000 photon packages per cell, one needs > 5 days of computational time on a modern 
PC (Pentium Xeon 2.4 GHz, 1Gb RAM). 

However, for certain cases the computational time needed for the ray tracing part can 
be efficiently reduced. To illustrate this for a rotating Keplerian disk with radial velocity 
gradient, let us consider a cell at a given radius r. The local line width that strongly depends 
on the temperature and microturbulent velocity in the cell is typically ~ 0.1 — 0.2 kms~^, 
which is significantly smaller than the total velocity dispersion across the disk, V^isp ^ 
10 kms~^. Therefore, due to the Doppler shift the photons emitted at a peak frequency in a 
cell rotating with a given velocity can only be absorbed locally by some of the neighboring 
cells having similar velocities (the so-called "interaction area"). An example of such an 
interaction area is shown in Fig. [2], where regions of the Keplerian disk that are radiatively 
coupled to the cell located at r 130 AU (white spot in the figure) are shown in black. 

We modified the ray-tracing algorithm used in the RT solver of "URAN(IA)" such that 
the iterative integration of the RT equation is only performed over the radiatively coupled 



http : //www . strw . leidenuniv . nl/ astrochem/ radtrans/ 
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regions of the disk. During the first iteration the radiatively coupled regions for each cell are 
identified using the full ray tracing and stored in a file. This information is utilized during the 
next iterations to integrate the RT equation over the radiatively coupled regions only. Note 
that the radiatively coupled region will be different for each line of each of the considered 
molecules. However, for those molecules for which the local lines widths are dominated by 
non-thermal motions (e.g., microturbulence) this difference will be small. 

This modification preserves the flexibility of the Accelerated Monte Carlo approach and 
allows to decrease the computation time for protoplanetary disks by about one order of 
magnitude and more, retaining the same high level of accuracy. One may use an analytical 
prescription instead of the preliminary ray tracing, but such a prescription will be model- 
dependent and thus cannot be easily generalized. Though our concept of the interaction areas 
in rotating configurations somewhat resembles the idea of the approximate LVG method (see 
Sect. 14. 3p . it leads to an exact solution of the LRT problem. 

3.2.2. Concept of thermalized cells 

The densities in the disk interior are so high that the low level populations of many 
molecules are thermalized. Naturally, the thermalized disk cells can be effectively excluded 
from the global iteration calculations since the radiative effects do not play a role in populat- 
ing the levels. We use the following approach to determine whether a disk cell is thermalized 
or not. First, the level populations in the disk are calculated using the simple LTE and Full 
Escape Probability (FEP) approaches (see next section). These two approximate methods 
represent two extreme cases to the solution of the LRT problem. The LTE approach tends 
to overestimate while FEP method tends to underestimate the populations of high levels in 
the case of radiative coupling. The cell is supposed to be thermalized if the relative differ- 
ence between the corresponding LTE and FEP level populations is smaller than an adopted 
threshold of 10^'^-10~^. This modification leads to an additional computational acceleration 
of about 10% for HCO+ and 30% for CO and CO isotopes. 

4. Approximate LRT methods 

In this section, we summarize the approximate LRT methods that will be compared 
with the results of the Monte Carlo simulations from the ART code. A schematic overview 
and comparison of the utilized LRT approaches is given in Table [31 
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4.1. Local thermodynamical equilibrium 



In the case of local thermodynamical equilibrium (LTE), the molecular level populations 
are thermalized. They follow a Boltzmann distribution with kinetic gas temperature Tkin'- 



Hi 



Hi 



ii 

9j 



exp 



kTi 



kin 



where gi and gj are the statistical weights and Uij is the frequency of the i — j transition. In 
terms of radiative transfer, thermalization means that the excitation and kinetic tempera- 
tures are equal for a given transition: 



kin 



T 



hu. 



ln(^^ 

9j n-i 



(2) 



When the excitation temperature is higher or lower than the kinetic temperature, the level 
populations are super- or subthermally excited. 

There are two general situations when LTE holds: 1) The hydrogen density is substan- 
tially higher than the critical density of a transition, or 2) the optical depth of a transition 
is so high that a perfect equilibrium between the collisional and radiative (de-)excitations is 
reached. The LTE a pproach is appropriate for the l ower r otational transition s of CO and 
its isotopologues (see Dartois. Dutrey. fc Guilloteau ( 2003 ): Pietu et al. ( 20051 )). It has also 
been used in the analysis of molecular line data even if the conditions for it s applicability 
are not fully iustified because it is the easies t LRT method to implement (e.g., 
2OO3I : bi et al.l[2003l : iNaravanan et al.lbood ). 



Aikawa et al. 



4.2. Full escape probability 



Another extreme is to assume that a transition has a very low optical depth, so that 
the newly generated photons escape the disk freely, without any further interaction with 
the medium. Or, by other words, in the PEP approach the escape probability /3 is not 
calculated and always set to 1 (see also Ivan der Tak et al.ll2007l ). Under such conditions the 
LRT problem becomes local and can be easily solved. The Pull Escape Probability (PEP) 
level populations are obtained as a solution of the statistical equilibrium equations with the 
mean intensity J substituted by the mean external radiation field, e.g. CMB: J = Jcmb- As 
in the LTE case, the PEP approach requires only very little computational time. PEP can 
be an adequate method for molecules with very low abundances. 
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4.3. Large velocity gradients 



In some objects such as expanding stellar envelopes or galaxies, radial veloci ty gradients 



can be much larger than the loc al thermal and microturbulent velocities (e.g., IWeiB et al. 



2005 



Castro-Carrizo et al.l 120071 ). In this case the photons that are emitted by a certain 
region can only be absorbed locally. The size of this region is determined by the local 
velocity field and the line width. If one assumes that the physical conditions in this local 
envi ronment are homogeneous, the radiative transfer prob lem can be substantially simplified 



le.g- 



Sobolevlll960l : ICastoi]ll970l : iGoldreich fc Kwanlll974[ ). 



The Keplerian rotation of a circumstellar disk causes strong velocity variations of ~ 
0.5-5 kms~^, which are larger than the local thermal (~ 0.2 kms~^) and microturbulent 
(~ 0.2 kms~^) velocities. As shown in Fig. [2|, the size of the radiatively coupled zone 
is significantly smaller than the entire disk plane. Therefore, one may expect the LVG 
approach to be efficient for the line radiative transfer modeling. 

In the LVG method, the mean intensity J in a disk cell is approximated by a combination 
of the source function S for a given transition and the external radiation field Jcxt'- 

J= (1-/5)5 + /3Jext, (3) 

where (3 is the probability for the photons to freely escape the medium, < (3 < 1. The 
escape probability f3 is usually expressed in the form: 

1 — exp(— r) 



P 



T 



(4) 



where r is the effective line optical depth because of the local velocity gradient. Note that 
for optically thin lines the escape probability /3 reaches unity and the LVG method turns 
into the FEP approach. 



The effective optical depth r is evaluated as: 

T = a X AL, 



(5) 



where a is the local absorption coefficient and AL is the coherence length i.e. the mean 
distance for a photon to escape the medium. 

The coherence length AL in the equatorial plane of the Keplerian disk can be roughly 
estimated as follows (see Fig. [3D. The photons emitted in a disk cell at a radius R cannot 
escape in the radial direction R since there is no velocity shift in this direction. Instead, 
we assume that the photons can only escape in the tangential direction. Ideally, one would 
need to apply spatial averaging to get a correct value of the coherence length, but this would 
significantly complicate the realization of the LVG method. 
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The approximate coherence length can be expressed via the local line width and the 
Keplerian velocity. The velocity shift along the coherence path is, by definition, equal to the 
local line width AV: 

Ay = |1^2Cos((^) - l^il, (6) 

where Vi and V2 are the Keplerian velocities at the distances of R and R + AR, respectively, 
and (p is the corresponding angle between the velocity vectors (Fig. [3]). Expanding the 
Keplerian velocity into Taylor series, we find that V{R + AR) V{R) + dV/dRAR, while 
cos{!f) ~ R/ {R + AR). This allows to derive the value of AR from Eq. ([6]): 

AR 2 AV 

IT = 3—- 

Note that AR is the width of the ring-like radiatively coupled region depicted in Fig. [2l 

From Fig. [3] we also find that the coherence length AL can be related to the radial shift 
AR and radius R as AR/ AL = AL/R, and thus 

AL = R^ll^ ^R'l' (8) 

for a constant line width AV . This approach is adopted in our LVG method. It typically 
takes only several seconds of CPU time to calculate the molecular level populations with the 
LVG approach. 



4.4. Vertical escape probability 



One has to keep in mind that in our disk models there is no velocity gradient in the verti- 
cal direction and the disk cells are radiatively coupled in vertical direction, while in the LVG 
approach, the disk is effectively decoupled in independent parallel planes. Let us compare 
the coherence length AL with the scale height of the disk Hj,. For the geometrically thin. 



isoth ermal disk Hp = \/2Cs/^ (for our definition of Hp, see IPartois. Dutrey. fc Guilloteau 



20031 ). where Cg is the sound speed and fl is the Keplerian angul ar velocity at rad ius R. For 
grey dust opacities the radial temperature profile is T oc (jDuUemondl I2OO2I ) and thus 

Hp oc -R^/^, similar to the scaling of the coherence length. For our disk model Hp ^ 28 AU 
at -R=100 AU, while AL ^ 25 AU. Therefore, the radiative effects in vertical and radial di- 
rections are comparable over the entire Keplerian disk (at least around the disk mid-plane) 
and, by design, the LVG method can only partly account for the radiative coupling. 



To better illustrate the role of radiative coupling in vertical and radial directions, we 
computed with ART the ratio between tangential and vertical optical depths for the uniform 
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disk model (see Fig. H]). The typical value of this ratio is only about a few over the entire 
disk, so both directions are important for the line radiative transfer. Using ART, we also 
calculated the mean coherence lengths in the equatorial plane of the disk. This is done 
by spatial averaging of individual coherence lengths. The calculated mean and individual 
coherence lengths are compared to the values from Eq. ([8]) in Fig. [51 As can be clearly seen, 
our analytical expression reproduces well the radial dependence of the mean coherence length, 
but a number of the individual lengths have larger values because of the photon trapping 
along directions with zero Doppler shifts. As a result, the mean value of the coherence length 
is typically about 10 times longer than predicted by Eq. ([8]). 

However, the excitation conditions vary so greatly over the mean coherence length that 
it cannot be directly put into Eqs. and (jl]). In order to properly evaluate the total 
escape probability (3 one has to average individual escape probabilities rather than the co- 
herence lengths. This requires information about exact physical parameters along all coupled 
directions, which makes this approach non-local and thus far too complicated. 

In the LVG approximation the coherence length is determined by the velocity gradients 
resulting from the Keplerian rotation of the disk. However, as we have shown above the 
radial radiative coupling is comparable to the coupling in vertical direction. One can assume 
that the vertical radiative coupling is more important than the radial coupling, which leads 
to significant simplification of the radiative transfer. 

This idea is utilized in the VEP method where photons are assumed to escape only in 
the vertical direction perpendicular to the disk plane. A simple correction is made to take 
into account the large asymmetry between the two opposite escape directions: the escape 
probability is computed by taking the sum of escape probability upwards and downwards: 

l/ l-exp(-r-) l-exp(-r+) ^^ 
^-2*^ + J' 

where r"*" is the optical depth towards the disk surface (upward opacity), and r~ is the optical 
depth through the disk midplane (downward opacity). The optical depth in upward and 
downward directions is calculated using the local excitation conditions and the molecular 
surface densities along the corresponding directions. Note that this assumption is rather 
coarse in the presence of strong vertical gradients of physical conditions. 

Similar to LVG, in the VEP approach the escape probability, source function, and 
the local optical depth are obtained iteratively. However, the direct implementation of the 
iterative VEP approach can be dangerous for disk cells with negative excitation temperature. 
In VEP these local conditions are extrapolated over the total vertical extent of the disk 
and the resulting optical depth can be negative and large. This leads to the physically 
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unrealistic situation that the escape probabihty can become large, especially during the 
iterations required to solve the statistical equilibrium equations and the VEP method may 
fail to converge. To avoid this problem, we control the downward line optical depth r~ such 
that it cannot become negative. The upward optical depth is allowed to remain negative, in 
order to represent weak inversions at the surface layers. 

The method was initially developed in the DiskFit software package (see section 4.6). 
Like FEP, the VEP method takes only a few seconds of CPU time to solve the LRT problem 
(see Table E]). 



4.5. Non-local ID-method VOR 

In the LTE and FEP methods radiative coupling effects are completely ignored, while 
in the LVG and VEP approach radiative trapping is only partially treated. Since there is no 
vertical velocity gradient in the adopted disk model, the cells are generally radiatively coupled 
in vertical direction. Furthermore, Fig I2] shows that, with the exception of a relatively narrow 
cone in the radial direction, a given cell in the disk is only coupled with cells at the same 
radius, i.e. with regions with the same ID structure. 

Using these ideas, we construct a non-local approximate LRT method in which the 
radiative coupling in only the vertical direction is taken into account ( "VOR" - Vertical One 
Ray method). This simplifying assumption allows to use a ID approach for calculating the 
mean line intensity J in a disk cell: 

^ = ^ (A + h) , (10) 

where Ii and I2 are the intensities of the radiation coming along vertical direction from the 
top and bottom of the disk, respectively. To further accelerate this method, only one central 
frequency for each emission line is taken into account. This non-local RT problem is finally 
solved using only 2 photon paths along vertical direction for the ray tracing. Computational 
demands of the ID VOR method are much less than those of the exact 2D code ART. It 
usually takes about a few minutes of CPU time to obtain the level populat ions. Note that 



our VOR approach is another ir nplementation of the Feautrier method (see lFeautrierlll964l ) 



and resembles the ID scheme of iLucad (119741 ). 



The ID methods have been used to represent spheres as well as infinite planes. In the 
latter case, an Eddington factor of 3 is used to take into account the mean opacity over all 
possible angles. In our case, since the interaction only occurs in a ring centered around the 
star, the Eddington factor should be smaller. We have simply used the Eddington factor of 
1, like for a sphere. 
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4.6. DiskFit 



The described approximate methods are included in the "URAN(IA)" package. In 
our study we also used another software package "DiskFit". The program DiskFit allows 
to model line radiative transfer in disks using various approximate methods as well as to 
reconstruct basic disk parameters using an advan ced x^-minimization scheme in the wf-plane 
(IDutrey et al.l Il994j : iGuilloteau fc Dutreyl Il998l ). The program offers three possibilities to 
solve the coupled equations of statistical equilibrium an d RT: LTE, VE P, and a non-local ID 



eautrier 



19641 ) and is an adaptation 



19741 ) for the line radiative transfer 



approach. The latter is based on the Feautrier s cheme (| 
to the disk geometry of the code developed by iLucad 
modeling of interstellar clouds. This ID method is thus similar to VOR, but uses several 
frequency points across the line profile. The Eddington factor may also be adjusted if needed. 

Once the level populations are derived, the interferometric visibilit ies as well as beara - 
convolved channel maps can be computed using the 3D ray tracing part (IDutrey et al.lll994l ). 
Any arbitrary orientations of the disk defined by the inclination and position angles can be 
used. For that, a three dimensional regular grid is used, with additional refinement for 
the central disk region where the surface brightness gradient is very strong. The DiskFit 
radiative transfer methods and the ray tracing part were thoroughly compared against those 
implemented in the "URAN(IA)" package and good overall agreement was achieved for the 
disk level populations and spectral maps for both the optically thin and thick molecular 
transitions. Therefore, for the rest of the paper we compare approximate and exact LRT 
methods from the "URAN(IA)" package only. 



5. Results 

Prior to the detailed comparison of various radiative transfer methods described above, 
we discuss the problem of molecular line formation and excitation conditions in disks in 
general. 



5.1. Molecular line formation in protoplanetary disks 

5.1.1. Excitation regimes with uniform abundances 
Since the th ermal and density structure of protoplanetary di s ks is in relatively we ll 



understood (e.g., iD'Alessio et al.l 1 19981 : iDuUemond fc TuroUal 12000 
we can discriminate distinct regimes of line formation. 



20011), 



D'Alessio et al 



- 16 - 



A particular rotational transition is excited when the molecular hydrogen density ex- 
ceeds a critical density ricr'- 

ncr = ^ • (11) 

Here A^i is the Einstein coefficient and Cut are the coUisional rates from the upper level u to 
all other lower levels i. The density, nth, at which the molecular levels become predominantly 
populated by collisions and thus are thermalized is typically about one order of magnitude 
higher than the critical density. We will further refer to nth as to the thermalization density. 

Consequently, the disk can be roughly divided in two regions. In the so-called "super- 
critical" region near the disk midplane, where the density is larger than the thermalization 
density, molecular level populations follow the Boltzmann distribution and the excitation 
temperature is equal to the kinetic temperature, Tex = ^l^in- In this region a simple LTE 
approach can be safely utilized. 

In the upper, more diffuse "sub-critical" region of the disk the line excitation is partly 
controlled by collisions and partly by internal/external radiation. The crucial parameter that 
determines the excitation conditions in this region is the molecular column density. If this 
column density is so high that the sub-critical region is opaque to the internal radiation, then 
the collisional and radiative (de-)excitations are in perfect balance and the corresponding 
molecular levels are thermalized, T^^ = Tkin. In contrast, when the column density is so low 
that the sub-critical region is fully transparent to the internal radiation, the corresponding 
molecular levels are mostly radiatively populated and therefore subthermally excited: Tex ~ 
Tbjf|. In this case a simple Full Escape Probability approximation can be used (see Section |5]). 
The characteristics of the emergent molecular line spectrum will depend on the relative 
contribution from each of the super-critical and sub-critical regions. 



5.1.2. Effects of chemical stratification 

This simple scheme becomes more complicated when the effects of chemical stratification 
are taken into account. In this situation three distinct regimes of line excitation can be 
distinguished. 

In the ffist, and most simple case the molecular layer is located deeply inside the super- 
critical density region and the molecular level populations are in LTE (see left panel in 
Fig. [6]). Good examples of that kind are the lower rotational transitions of CO and its 



■^Here Tbg is the temperature of the background radiation, e.g. 2.73 K for the CMB field 
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isotopologues (jPartois. Dutrey. fc Guilloteaij|2003l ). 



In the second case, the molecular layer is entirely confined within the sub-critical density 
region. The distribution of the level populations will then depend on the molecular column 
density as discussed above (middle panel, Fig. [6]). When this column density is low and 
radiative coupling is not important, the FEP approximation works. In the opposite case the 
radiative transfer problem has to be tackled with a more sophisticated no n-local method. 



The line excitation proceeds in this way in gas-deficient transition disks (jjonkheid et al. 



20061). 



Finally, in the most general case the molecular layer is located between the sub-critical 
and super-critical disk regions (right panel in Fig. [6]). In this situation the LTE approach does 
not work, whereas the FEP approximation is valid only if the optical depth of the sub-critical 
layer is low and the molecular emission is dominated by the super-critical (thermalized) layer. 
The lower rotational lines of e.g. HCO^, HON, and CS are likely excited in such a way. 

If the molecular column densities are so large that the sub-critical region becomes op- 
tically thick for internal radiation, the simple FEP approach fails and non-local radiative 
transfer modeling is required. The non-LTE effects in disks are likely important for the 
rotational lines of complex species, e.g. H2CO, but may also be important for higher tran- 
sitions of many other simple molecules (such as HCN, CS, etc.). In the next section we 
verify whether non-LTE effects are important in protoplanetary disks and for which molecu- 
lar transitions this is the case. There we use a simple preliminary analysis of the excitation 
conditions. 



5.2. Critical column density diagrams 

It is often useful to understand under what conditions a given transition is excited be- 
fore any LRT modeling is performed. One example is the iterative x^-minimization analysis 
of interferometric data where an approximate LRT method has to be utilized due to lim- 



ited c omputational resources (e.g.. lGuilloteau fc Dutreylll998l : lDartois. Dutrey. fc Guilloteau 



20031 ). Usually the basic properties of the physical structure and sometimes the chemical 
structure of the disks are at least partly known. When the chemical structure is not known, 
one can use uniform abundance distributions with the abundance values taken from theo- 
retical models. 

Using the density and chemical structure of the disk and molecular data as input, one 
can determine the line optical depth and the amounts of emitting molecules in the sub- and 
super-critical regions. As discussed in Section 14. 5^ radiative coupling is strong in vertical 
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direction due to the absence of velocity gradients, and thus radiative effects on the hne 
excitation are most profound for face-on oriented disks. Plotted as a function of radial 
distance, the molecular column densities in the sub- and super-critical disk regions as well 
as the molecular column density needed to reach an optical depth of unity may serve as 
an indicator whether and where the considered transition is non-LTE excited and which 
approximate or exact method is appropriate to model the radiative transfer in this case (see 
previous section). 

In Fig. [7]we present such Critical Column Density (CCD) diagrams for some transitions 
of the molecules under investigation: CO, CS, HCN, H2CO, and HCO"^. We focus on low 
rotational transitions and adopt approximate values of their thermalization densities, viz., 
10^ cm~^ for the CO lines and 5 x 10^ cm~^ for all other molecules. The column densities to 
attain r = 1 are computed for the central velocity channel of the (1-0) transition (Iqi — Oqo 
for para-H2C0) assuming a constant excitation temperature of 10 K. 

Not surprisingly, the column density of CO is so high, A^(CO) ~ lO^^cm"^, that the 
lines are very optically thick, r ~ 10^-10^. Since the CO lines are easily excited already at 
densities of about 10^ cm~^, essentially all CO molecules are located within the super-critical 
region (case "a)" in Fig. [6l) . Molecules tracing the disk interior are excit e d in similar nianner , 
e.g. H2D+, N2H+, DCO+ jAikawa fc HerbstlboOll : Isemenov et allbood : butrev et ahlboO?! ). 
Note that in our model, the CO abundances are high (see Table 2). When comparing with 
real disks, where CO can be significantly depleted, our test case for C^^O can be applied 
to ^^CO observations, while for ^^CO, the LTE approximation will work even better (since 
r = 1 will be reached at higher densities). 

For other, less abundant species that have higher critical densities and are concentrated 
in the disk intermediate layer, the CCD diagrams show a bimodal behavior. In the inner, 
dense disk region at r < 400-600 AU almost all molecules are located in the super-critical 
region, and their levels are LTE-populated (similar to CO). In the outer, less dense disk 
region the molecules are mainly concentrated in the sub-critical layer (case "b)" in Fig. [6]). 
The optical thickness of the sub-critical region for the disk model with chemical stratification 
is lower or close to one. Consequently, the molecular sub-critical layers are expected to be 
optically thin or moderately optically thick in the outer disk region and the FEP approxima- 
tion holds (see HCN and HCO"'" in Fig. [7]). Thus, non-LTE effects should not play a major 
role for low rotational lines of most molecules in protoplanetary disks. A notable exception 
is water that has a comple x level structure wit h a number of maser tra nsitions and some 
highly optically thick lines (jPoelman et al.l 120071 : Ivan der Tak et al.l |2006| ) . 



In contrast to the layered chemical structure, for the uniform model the amount of 
molecules is so high in the sub-critical region that this region is optically thick and radiative 
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trapping becomes efficient. As a representative example, we present CCD diagrams for the 
uniform HCO^ model (lower right panel in Fig. [7]). In what follows, we will focus on this 
uniform HCO"'" model to analyze in detail the importance of non-LTE excitation for both 
low and high transitions as modeled by various LRT methods. 



5.3. Excitation temperatures 

In Fig. [8] we show the spatial distributions of the excitation temperature in the uniform 
disk model as computed by the different LRT methods for the HCO+ (1-0), (4-3), and (7-6) 
transitions. By definition, the LTE excitation temperature for any transition is equal to the 
gas kinetic temperature. The vertical temperature gradient and cold isothermal midplane 
are clearly visible (Fig. [HI top row). In contrast, the excitation temperature distributions 
calculated with other LRT methods show a departure from LTE conditions, in particular for 
high transitions and for the upper, less dense disk regions. 

The most extreme deviation from LTE is the negative excitation temperature caused 
by the inversion in the populations of the first two rotational levels of HCO^ (maser effect), 
which is reached in the warm upper disk region, at r < 200 AU (white zone in Fig. [8|, left 
panels). In this part the (self-)radiation is not intense, and the levels are excited and de- 
excited by collisions but only radiatively de-excited. Consequently, LTE is not longer valid 
and coUisional pumping of the rotational level populations appears as a result of a specific 
ratio between radiative and collisional (de)-excitation probabilities. Such inversion is char- 
acteristic of some hnear mol ecules like CS, HCO"* " , and SiO and requires high temp eratures 



and a specific density range (ILiszt fc LeungI 119771 : IVarshalovich fc Khersonskiilll978l ). 



However, this inversion does not result in maser amplification of the HCO"'" (1-0) line 
intensity even for highly inclined disks. As soon as the line opacity and thus the self-radiation 
becomes significant, the stimulated radiative transitions depopulate the upper levels and 
destroy the inversion. Note that in the FEP approximation the effect of self-radiation is 
completely neglected and therefore the negative excitation temperature zone is strongly 
overestimated (Fig. [HI left bottom panel). 

Another non-LTE zone corresponds to the subthermally excited sub-critical disk region, 
which occupies a significant disk volume for high HCO"*" transitions. 

To quantify the difference between the excitation temperature T^x computed with an 
approximate and the exact ART methods, we introduce the local error e in each disk cell: 

T — T^f^-T 

(12) 



|7^ex| + |TART|' 
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A relative difference of 20% in the excitation temperature corresponds to e ~ 0.1, and a 35% 
difference leads to e ^ 0.2. 

Consequently, the global accuracy Y of an approximate method over the entire disk can 
be evaluated with the following expression: 



where the absolute values of the local errors ej are summed over the optically thin disk cells 
and weighted by the number of molecules Mj in each cell. M is a total number of molecules 
in the optically thin region. The optical depth at the line central frequency, tj, is calculated 
in the vertical direction using the ART level populations. 

This global accuracy criterion is designed to directly relate the difference in the excita- 
tion temperatures to the difference that will appear in the spectra. We consider disk cells 
with r < 1 (in vertical direction) because this is the region where most of the line emission 
comes from (for face-on disks). In essence, the errors in the r ^ 1 regions are not important, 
since these regions are not visible. Certainly, this criterion is not a unique one and has to 
be used with care. 

We state that an approximate method gives "good" accuracy ( "A" ) if the global criterion 
Y < 0.1, "reasonable" accuracy ("B") when 0.1 < F < 0.2, and "bad" accuracy ("C") 
when Y > 0.2. The distributions of the local error e and global accuracy Y of the applied 
approximate LRT methods for the uniform HCO"'" model are shown in Fig. [91 

All LRT methods are accurate in the super-critical disk region, n > 10^ cm~^. However, 
the thermalized disk region extends toward lower densities because the HCO"'" lines are highly 
optically thick and the perfect balance between radiative and coUisional (de-)populations is 
reached (see Sect. 15.1.11) . The thermalized region is particularly large for low rotational 
transitions that are excited at lower densities. 

The LTE approximation tends to significantly overestimate the excitation temperatures 
for the high transitions with high critical densities but is still appropriate for the J =(1-0) line 
(Y < 0.1, Fig. [9]). In contrast, FEP fails for the (1-0) transition since it greatly overestimates 
the region of negative excitation temperature. It is rather accurate {Y ^ 0.16) for the (4-3) 
line because it adequately describes the excitation conditions in the small, optically thin 
zone above r = 1 in the outer disk region, at r > 400 AU, which occupies the largest volume 
in the disk. FEP is also accurate for (7-6) transition {Y = 0.12), since it works in the large 
optically thin region that is decoupled from the internal radiation. 

The LVG, VEP and VOR methods are in general much more accurate than LTE and 




(13) 
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FEP because they account for the radiative trapping. In particular, the maser zone in the disk 
is well reproduced by LVG and VOR. The LVG approach shows the same tendency for the 
excitation temperature deviations as FEP because it underestimates the radiative coupling. 
The non-local VOR approach gives the most accurate results for all considered transitions 
since it treats the radiative transfer in vertical direction, where velocity gradients are absent 
in our model. A minor problem of VOR is the overestimation of the photon trapping in the 
very upper disk regions where the radiation can escape also in other directions. 

In general, for the uniform HCO^ model and all considered transitions FEP and LTE are 
not accurate LRT methods ("C" grade), while LVG and VEP are acceptable ("B"), and VOR 
is good ("A"). This conclusion is valid in general for all other key molecules apart from CO 
isotopologues where all methods are accurate because their lines arise from the super-critical 
region (see Table HI). For the disk model with chemical stratification the overall accuracy of 
all LRT methods is better because the lines are more optically thin and radiative effects are 
not that important (Table [5]). For many molecular transitions LTE and FEP are reliable 
approximate methods for the line radiative transfer modeling of chemically stratified disks. 
However, in this case both LVG and VEP provide an even better approximation, at the 
expense of very little extra computational time. 



5.4. Synthetic spectra and spectral maps 

Using the excitation temperatures calculated for the uniform chemical structure, we plot 
the corresponding synthetic HCO"*" spectra for the J = (1 — 0), (4-3), and (7-6) transitions 
in Fig. [TDl All spectra are convolved with the same 10" beam and are computed for two disk 
inclinations: i = 0° (face-on) and i = 60°). 

The LVG, VEP and VOR methods result in line profiles with intensities deviating by 
no more than 30% from ART, whereas the FEP and LTE intensities and line profiles differ 
more strongly from those of ART, in agreement with the difference in excitation tempera- 
ture maps (see Fig. [HI and Tables HIS]) • The line intensities computed with FEP and LTE 
are either overestimated or underestimated, depending on transition. The overall disagree- 
ment between these approximate methods and ART does not depend strongly on the disk 
orientation. 

For the HCO+(1-0) transition the LTE spectrum is close to the ART line profile, but the 
FEP intensity is 3 times higher for the face-on orientation and even higher for the inclined 
disk. As discussed in the previous section, this is due to the maser effect in the first two 
HCO"'" level populations. The same effect applies to the lower transitions of CS (not shown 
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here). However, for the higher HCO"'" transitions the FEP spectra are in better agreement 
with the full Monte Carlo simulations, whereas LTE tends to significantly overestimate the 
line intensities. This is because the size of the radiatively coupled part is decreasing for 
higher transitions, having higher thermalization densities and thus a more extended sub- 
critical region where Tex ^ T\^in (compare top panels in Fig. [9]). 

At face-on disk orientation a self-absorption dip arises in the (4-3) and (7-6) ART spectra 
at the central velocity channels. This dip is caused by absorption of emitted radiation in 
the upper disk region, at n < lO^cm"^, due to the vertical excitation temperature gradient 
(Fig. [8]). Both LTE and FEP fail to accurately reproduce this effect. In contrast, the 
double-peaked HCO+ spectra of the inclined disk are not due to self-absorption but appear 
as a result of the beam convolution over disk cells with different projected velocity. In this 
case all approximate methods correctly predict the shape of the line profiles. Note that while 
integrated line intensities do not strongly depend on the disk orientation, the peak intensities 
do differ strongly, and are lower for higher transitions. 

The effects of chemical stratification on the integrated spectra are shown in Fig. [TTl 
Since HCO^ occupies only a fraction of the entire disk, the computed line intensities and 
their optical depths are much lower compared to the uniform model. The self-absorption 
dip is absent in the face-on line profiles since there is no more HCO"'" in the upper disk 
region. The overall agreement between the different LRT methods becomes better (see also 
Table E]). Still, LTE fails for higher transitions. The detailed comparison for the other species 
considered here yields similar results. 

The criterion Y that gives a global accuracy estimate for the excitation temperature 
calculated with an approximate LRT method is also valid for the single-dish spectra. To 
investigate how good this criterion is for the synthetic spectral maps, we plot the HC0~*"(4- 
3) map at the 0.68 km s~^ velocity channel for the same 60°-inclined disk, but without beam 
convolution (Fig. [T2l). 

The model with the uniform HCO^ abundances leads to a fan-like structure with two 
low-intensity lanes and a high-intensity interior zone (see Fig. [T21 upper row, central panel). 
The low-intensity lanes correspond to the cold disk midplane with low excitation temper- 
atures. Similar to the single-dish spectrum (Fig. [10]), LTE (FEP) tends to overestimate 
(underestimate) the total map intensity. The LTE significantly overestimates the excitation 
temperature in the upper disk region, leading to the formation of two high-intensity features 
near the vertical edges of the image. The contrast of the FEP spectral image is much smaller, 
yet the cold midplane and high-temperature inner zone are clearly visible in the map. 

The LTE, FEP, and ART HCO^(4-3) maps simulated with the layered chemical struc- 
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ture do not differ so much as for tlie uniform model (see also Fig. [TT]) . However, the spectral 
images show a more complicated pattern that resemble two thick elliptical rings overlaid on 
each other with an angular shift of about 20° and fixed at the center. The three low-intensity 
lanes correspond to the disk region close to the midplane that is devoid of HCO"*" emission 
at the chosen velocity channel of 0.68 kms~^. Furthermore, the HCO"'" ions are also absent 
in the disk upper layer and the overall extension of the map is smaller in vertical direction 
compared to the uniform model. 



6. Discussion 

Our analysis of the hne r adiative transfer in protop lanetary disks supports and ex- 



tends the analysis performed by Ivan Zadelhoff et al.l (l200ll ). In particular, they argued that 



SE( Ji,=0) calculations (in our notation FEP) provide better agreement with SE calculations 
(in our notation ART) compared to LTE. According to our results, the FEP method is 
ind eed more correct th a n LT E for the upper molecular transitions which were considered 



by Ivan Zadelhoff et al.l (1200 ll ). However, the FEP approach may give completely wrong 
intensities for the lowest molecular transitions such as J = 1 — and J = 2 — 1 since it over- 
estimates the maser effect, resulting from the specific ratio between coUisional (de)excitation 
rates. Hence, the LTE approach provides more reliable intensities of the lowest transitions 
than the FEP method. However, LTE fails to reproduce important features of the line pro- 
files like the self- absorption dips. Despite the fact that either LTE or FEP can be sufficient 
for the simulations of the particular molecular lines, we recommend to use LVG or VEP 
since they provide in general better agreement at the same cost of complexity and CPU 
requirement. The more exact, but more time-consuming methods such as VOR or ART can 
be used to check the validity of the above approximate methods. 

An interesting result of our study is the fact that the synthetic images calculated with 
VEP are quite close to the exact ART spectra, for both chemically uniform and layered 
models. This may look surprising since the idea of the VEP method is to use local exci- 
tation conditions and to extrapolate them along the vertical direction to calculate escape 
probabilities. Thus, one may expect that this "homogeneous" approach is inappropriate for 
the disk models with their strong gradients of physical conditions and molecular abundances. 

The "good behavior" of VEP is a consequence of the very strong (vertical) density 
gradient compared to the temperature gradient. In the layered disk model molecules are 
concentrated in thin layers where physical conditions do not change much in vertical direc- 
tion. In the uniform disk model the strong vertical density gradients result in a super-critical 
region around the midplane where molecular emission is thermalized, while emission from 



-24- 



a low-density surface layer is negligible. In between lies a narrow sub-critical disk region, 
where molecular opacities can be high enough to affect the emergent spectra (see Sect. 15. ip . 
For optically thin lines, the contribution of this narrow layer to the emission is small. For 
optically thick lines, the likelihood that the line opacity reaches about unity in this region is 
very low. The exact location of the r = 1 layer determines the brightness temperature, but 
is not very critical since vertical temperature gradients are not very large. 



Recently iPoelman fc Spaand (120051 ) have proposed a new approximate Mu lti-Zone Es- 

cape Probability (MEP) method to model water line emission from a PDR region (IPoelman fc Spaans 
20061 ). In contrast to our VEP approach, in MEP the probability for photons to escape a 
cell is computed in 3D, by spatial averaging of the optical depths over several (6) directions. 
Due to its 3D nature, MEP may treat more consistently the photon escape from complicated 
radiatively coupled regions in Keplerian disks (see Fig. [2]) than our Vertical Escape Proba- 
bility method. Since the ID VEP approach already reasonably accurate in most cases, the 
use of 3D MEP for the LRT modeling of protoplanetary disks is recommended in situations 
when even higher accuracy is required at the expense of extra CPU costs for the spatial 
averaging. 



Another recent ID escape probability method has been derived by lElitzur Sz Asensio Ramos 



( 120061 ). Their Coupled Escape Probability (CEP) approach is based on mathematical trans- 
formation of the coupled system of the integro-differential RT equation and linear balance 
equations into one system of implicit non-linear equations. Due to fast convergence of this 
implicit scheme, CEP gains about 1-2 orders of magnitude in computational speed compared 
to conventional ALI schemes. However, this method cannot be easily adapted to multidi- 
mensional LRT problems with multi-level molecules because of the mathematical "trick" 
that lies behind its algorithm. Taken at face value, in ID geometry, it would be conceptually 
similar to, but more sophisticated than the ID Feautrier approach implemented in our VOR 
method. 

In the present study we did not consider molecules with complex level structure, e.g. 
molecules with hyperfine com ponents since they require more detailed investigations. For 



instance, iDaniel et al.l (120061 ) used new collisional rates to simulate the N2H"'" lines from 
prestellar cores and showed that the usual assumption that N2H"'" sub-level populations 
follow the Boltz mann distribution may l ead to significant errors in the intensities of the 
hyperfine lines. Ivan der Tak et al.l (120061 ) mentioned severe problems in LRT simulations 
of II2O lines. Apparently, the formation of such complex molecular lines in protoplanetary 
disks requires a more extended analysis which has to be based on the appropriate collisional 
rates. 



In the present paper, we also did not treat the effects of dust emission on the line 
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formation. For the protoplanetary disks the dust emission at milhmeter wavelengths is 
neghgible compared to the molecular transitions. However, at shorter wavelengths the dust 
emission may be sufficiently strong to excite upper molecular transitions and, therefore, 
should be considered. 

Nor did we include the effects of vertical mixing and turbulence in our model. From a 
chemical point of view, turbulent diffusion acts as an efficie nt non-thermal desorption mecha 



nism and smoothes molecular abundance gradients in disks (jllgner et al.ll2004j : ISemenov et al. 



20061 : IWillacy et al.ll2006l : iTscharnuter fc Gaill 120071 ). Consequently, the results for such a 
disk model would be in between the results for our uniform and layered chemical models. 



7. Conclusions 



We analyzed the line radiative transfer and the formation of rotational lines of several 
key molecules (CO, CS, HCO"*", H2CO, HCN) in protoplanetary disks. A detailed model of 
the disk physical structure with vertical temperature gradient and uniform/layered chemical 
structure is used. We used the exact Accelerated Monte Carlo code ART and compare the 
results with several approximate LRT approaches, including LTE, the Full Escape Probability 
(FEP), LVG, and the Vertical Escape Probability method (VEP) and a ID non-local code 
VOR. Our main conclusions are as follows: 



1. Prior to any modeling of the radiative transfer, one can roughly estimate the im- 
portance of non-LTE effects for the line excitation and thus choose an appropriate 
approximate or exact method by using Critical Column Density (CCD) diagrams. The 
basic idea of the CCD diagrams is to calculate and compare the amount of emitting 
molecules in the sub- and super-critical regions, taking into account the optical depth 
of the line. With CCD diagrams we found that the upper disk layers are subthermally 
excited for many molecular lines, in particular for the upper transitions. 

2. By comparing excitation temperatures, synthesized spectra, and spectral maps we 
demonstrate that the simple LTE and FEP methods work reasonably well for chemi- 
cally stratified disks where molecules are located mainly in the warm intermediate layer 
or in the disk midplane, but are not always accurate for disk models with uniform abun- 
dances. In contrast, the LVG, VEP, and especially VOR methods take the radiative 
coupling in the disks into account and thus accurately reproduce the line intensities 
and profiles. Since the LVG and VEP methods are computationally inexpensive, their 
use is recommended. 
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3. We found a possibility to accelerate the ray-tracing part of the RT solver in Monte- 
Carlo method for rotating disks. In this modification the ray tracing is only performed 
over radiatively coupled disk zones, leading to a computational speed gain by factors 
of 10-50. 

This work is only the first step in our systematic analysis of the line radiative transfer 
in protoplanetary disks. The excitation of molecules with complex level structure, including 
hyperfine line components, and the applicability of various approximate LRT methods for 
the x^-minimization fitting of the observational data will be addressed in future studies. 

We thank Michiel Hogerheijde and Jinhua He for useful and stimulating discussions. 
Authors are thankful to the anonymous referee for valuable comments and suggestions. This 
research has made use of NASA's Astrophysics Data System. 
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Fig. 1. — Left: The distributions of density, temperature, and radial velocity in the adopted 
disk model. Right: The distributions of the molecular abundances relative to the total 
amount of hydrogen nuclei for CO, HCO"*", and CS. 
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X/R, 

Fig. 2. — The equatorial plane of a Keplerian disk with radius Rq. The white spot represents 
a rotating cell that emits radiation at a certain rest frequency. The regions where this radi- 
ation can be absorbed by the same molecules are shown in black (the so-called "interaction 
area"). The interaction area of the disk has velocities that are Doppler-shifted from the 
velocity of the emitting cell by no more than the local line width (expressed in kms~^). The 
ring bounded by white circles corresponds to the radiatively coupled region in LVG method. 
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Fig. 3. — Analytical estimation of the coherence length for a rotating disk. 
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Fig. 4. — Ratio of tangential to vertical optical depths calculated with ART for the uniform 
disk model and the HCO+(4-3) hne. 
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Fig. 5. — Radial dependence of the coherence length (Eq. [HI thick solid line) compared to 
the individual (circles) and 3D spatially-averaged (thin solid line) coherence lengths (ART). 
The uniform model and the HCO"''(4-3) line are used. 
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Fig. 6. — Distinct cases of line formation in protoplanetary disks with chemical stratification. 
Left: The molecular layer is located inside the super-critical (high) density region adjacent 
to the midplane and thus LTE holds for the level populations. Middle: The molecular layer 
is located within the sub-critical (low) density region close to the disk surface. The level 
populations can be accurately calculated either with the FEP method when the molecular 
column density is sufficiently low or require a more sophisticated non-local LRT method 
(e.g. VOR) when the molecular column density is so high that radiative coupling within the 
disk becomes important. Right: In the most general case the molecular layer lies both in 
the sub- and super-critical regions. Consequently, FEP is only accurate when the emission 
mostly comes from the super-critical region of the molecular layer and the optical depth of 
the sub-critical part is low. In all other situations non-local line radiative transfer methods 
have to be applied. 
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Fig. 7. — (From top to bottom, from left to right) Critical Column Density (CCD) diagrams 
for CO, CS, HCN, H2CO, and HCO^. For each of these species, the CCD diagram shows the 
radial distribution of 3 column densities in the chemically stratified disk: the column density 
in the sub-critical layer (solid line), in the super-critical layer (dashed line), and the column 
density that is needed to reach an optical depth of unity for the (1-0) transition (Iqi — Oqo 
for H2CO) (r ~ 1, thin solid line). The uniform HCO"'' model is shown in the lower right 
panel. Note that for better presentation the column densities are restricted to values above 
PS 10^ cm~^. The ratio between the column densities in the sub- and super-critical layers, 
together with the value for r = 1 , can be used for preliminary analysis of the line excitation 
conditions in disks and stimulate the choice of a proper LRT method. 
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Fig. 8. — (From left to right) Distributions of excitation temperature (gray scale, in Kelvin) 
for the HCO+ J =(1-0), (4-3), and (7-6) transitions computed with the uniform disk chem- 
ical structure (polar coordinates). The black line depicts the location of r = 1 for the 
corresponding transition, as seen from above the disk in vertical direction. Vertical scale is 
angle theta in radians. Shown are the results for the LTE, VOR, ART, LVG, VEP and FEP 
radiative transfer methods (from top to bottom). The white contours correspond to number 
densities of 10^, 10^, and 10^ cm~^. The white area in the (1-0) distributions marks the zone 
of negative excitation temperatures (maser effect). 
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R. AU 

Fig. 9. — The same as in Fig. [8] but for the local accuracy estimate e. The thick black line 
depicts the location of r = 1 for the corresponding transition, as seen from above the disk 
in vertical direction. The accuracy of an approximate LRT method is good when the global 
error Y < 0.1, moderate when 0.1 < F < 0.2, and bad when Y > 0.2. 
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Fig. 10. — (From top to bottom) Synthetic spectra for the HCO"^ (l-O); (4-3), and (7-6) 
transitions convolved with a 10" beam. The distance to the disk is 140 pc. The results 
are obtained with the 2D exact ART (thick solid line), the ID non-LTE VOR (dotted line), 
LTE (dashed hne), LVG (tiny solid line), VEP (dash-dotted hne) and FEP (thin solid hne) 
approaches using the disk model with uniform HCO''" abundances. The disk inclination is 
0° (left panels) and 60° (right panels). The intensity is given in units of the main beam 
temperature (Kelvin). 
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Fig. 11.— 
model. 



The same as in Fig. [TO] but for the layered HCO^ abundances from the chemical 
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Fig. 12. — HCO"'"(4-3) intensity maps synthesized without beam convolution for the 
0.68 kms^^ velocity channel and for a disk inchnation of 60°. The distance is 140 pc. 
The results are simulated using the level populations computed with the LTE (left), ART 
(middle), and FEP (right panel) methods. The uniform (top row) and layered (bottom row) 
abundances of HCO"*" are utilized. The intensity is given in units of the radiative temperature 
(Kelvin). 
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Table 1. Adopted parameters of the disk and central star 



Parameter 


Symbol 


Value 


Distance 




140 pc 


Stellar temperature 




4 000 K 


Stellar radius 




2MRq 


Stellar mass 


M, 


0.7 Mq 


Disk inner radius 




0.037 AU 


Disk outer radius 


-Rout 


800 AU 


Disk mass 


Mdisk 


0.07 Mo 


Surface density at 1 AU 


Si 


pa 63* gcm~^ 


Density profile 


p 


~ 1 * 



*Density structure is outcome of disk model. 
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Table 2. Parameters of the disk chemical model 



Molecule 




Molecular layer 


Abundance, 


Abundance, 




lower bound 


upper 


bound 


iV(X)/A^(H) 




^min 


^min 




^max 


(layered) 


(uniform) 


CO 


0.7 


0.0 


2.0 


0.0 


8 (-05)" 


1 (-04) 


Ci^O 


0.7 


0.0 


2.0 


0.0 


2 (-07) 


2 (-07) 


HCO+ 


1.37 


-0.1 


1.68 


-0.05 


6 (-10) 


1 (-08) 


DCO+ 


0.0 


0.0 


0.4 


0.0 


5(-ll)'' 


2 (-10) 


H2CO 


1.37 


-0.1 


1.2 


0.0 


3 (-10) 


1 (-08) 


HCN 


1.75 


-0.16 


1.92 


-0.07 


1 (-09) 


1 (-09) 


CS 


1.44 


-0.07 


1.83 


-0.03 

b 


1 (-09) 


1 (-08) 



^The value A{-B) reads as A x 10"-^. 

is located in cold outer disk regions at r > 100 AU. 
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Table 3. Overview of the applied LRT approaches 



Method 


Type 


AppUcability 


CPU time* 


LTE 


Local 


Super-critical 


< 1 s 






region 




FEP 


Local 


As above & optically 


< 1 s 






thin lines 




LVG 


Local 


Regions with strong 


~ 1 s 






velocity gradients 




VEP 


Local 


Regions with low 


~ 1 s 






velocity gradients 




VOR 


Non-local 


Regions with low 


15 m 






velocity gradients 




ART 


Exact 


Everywhere in disks 


24 h 



*CPU times are given for the 128 x 128 disk grid, 200 
photons, and a Pentium 4 2.4 GHz PC. 
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Table 4. Reliability of the approximate LRT methods for uniform model 



Molecule Methods 





LTE 


VOR 


LVG 


VEP 


FEP 


CO 


AAB* 


AAA 


AAA 


AAA 


ABB 


cs 


ACC 


AAA 


ABA 


ABA 


BCB 


HCO+ 


ACC 


AAA 


BAA 


CBA 


CCB 


HCN 


ACC 


AAA 


ABA 


BAA 


CCB 


H2CO 


ACC 


AAA 


AAA 


AAA 


BCA 


C^^O 


AAA 


AAA 


AAA 


AAA 


AAA 


DCO+ 


ABC 


AAA 


AAA 


AAA 


AAA 



*"A" - good accuracy {Y < 0.1); "B" - moderate 
accuracy (0.1 <Y < 0.2); "C" - bad accuracy {¥ > 
0.2). For all species except H2CO the 1st, 2nd, and 
3rd letters correspond to the (1-0), (4-3), and (7- 
6) transition, respectively. For H2CO the loi-Ooo, 
4o4-3o3, and 423-322 para-transition are considered. 
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Table 5. Reliability of the approximate LRT methods for chemical model 



Molecule Methods 





LTE 


VOR 


LVG 


VEP 


PEP 


CO 


AAA* 


AAA 


AAA 


AAA 


AAA 


cs 


AAC 


AAA 


AAA 


AAA 


AAA 


HCO+ 


ABC 


AAA 


AAA 


AAA 


BBA 


HCN 


ACC 


AAA 


AAA 


AAA 


CBA 


H2CO 


BBC 


AAA 


AAA 


AAA 


AAA 


Ci^O 


AAA 


AAA 


AAA 


AAA 


AAA 


DCO+ 


AAA 


AAA 


AAA 


AAA 


AAA 



*See notation in Table |H 



